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ABSTRACT 



A computer program was coded to obtain static aeroelastic effects 
on simple planforms through the use of subsonic lifting surface theory. 
The program was divided into two major technical areas, aerodynamic and 
structural, with matrix notation used to indicate the influence of each in 
the final aerodynamic loading distribution. Several typical stability 
derivatives were then obtained to make an accurate stability and control 
analysis on the desired planform. 
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I . INTRODUCTION 



In making a stability and control analysis on modern high perform- 
ance aircraft, the ability to make span load evaluations on the wing and 
tail with static aeroelastic effects taken into account is most important. 
It seemed desirable to develop a simple computer program utilizing 
current 'state-of-the-art' theories, and proper documentation, so that 
design groups might have readily usable methods at hand. 

Span load determinations have gone through a long line of 
evolution, right up to the lifting surface theories of today. Prandtl's 
lifting line theory is considered as the earliest, which was suitable 
for simple planforms of high aspect ratio and unswept quarter chords 
operating in incompressible flow. Modifications of this method were 
made by Glauert, Multhopp, Weis singer, and others until the present 
lifting surface theories suitable for low (or high) aspect ratio plan- 
forms with swept quarter chords operating in compressible flow were 
attained . 

Two lifting surface concepts for subsonic calculations are: 

1. Define the surface pressure loading in the chordwise and span- 
wise directions by admissible sets of orthogonal functions with 
the magnitudes of the unknown coefficients being obtained by 
collocation techniques. Ref. 6. 

2. Define a grid of boxes over the wing, each of which contains 
a horseshoe vortex, Ref. 4. Then determine the strengths 
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of the vortices by imposing the constraint of flow tangency at 
an equal number of control points . 

The approach chosen for this study was the second, which is 
usually denoted as a vortex lattice method and is thoroughly described 
by Belotserkovskii , Ref. 4. The major advantage to this approach 
is the convenient use of matrix methods for taking the static aero- 
elastic effects into account. 

The method of Gray and Schenk, Ref. 1, used initially on the 
Boeing B-52 aircraft, was carefully reviewed as background for this 
study. Their analysis, which was restricted to aspect ratio wings of 
greater than five, employed the Weissinger L-Method to account for 
subsonic aerodynamics. The advantage of the lifting surface theory 
used for this study is that wings of lower aspect ratio may be ana- 
lyzed and the ability to more accurately predict chordwise (as well 
as spanwise) distribution of aerodynamic loading is provided. The 
method of Ref. 1 did not allow for aerodynamic loadings other than 
additional type/ whereas the method of this study is not restricted 
to loadings with the aerodynamic centers at the local quarter chord. 

To account for compressibility in subsonic flow the Prandtl- 
Glauert rule was used, Ref. 2, sometimes referred to as planform 



Additional type pressure loadings are defined as the air loadings 
due to angle of attack only. Generally they have a cotangent type chord- 
wise pressure distribution in subsonic flow with a logarithmic type singu- 
larity at the leading edge, and vanish at the trailing edge to satisfy the 
Kutta condition. 
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distortion. The aeroelastic effects will depend upon the original wing 



structure and will vary proportionally with the free stream dynamic 
pressure. The span load evaluation at zero dynamic pressure corres- 
ponds to the rigid wing results. Symmetric span loads will be appli- 
cable to longitudinal stability and control analyses, with anti- symmetric 
to lateral-directional. 

The program has been restricted to simple planforms with 
straight taper, and the influence of fuselage interference was not 
taken into account. This was not an oversight but deliberate since 
it was felt appropriate at this time to obtain a straight-forward docu- 
mentation on a simple wing situation, with applications to more com- 
plex planforms and interference effects as a logical extension of this 
study. 

The approach used in the structural area was restricted to a 
basic beam-like wing structure, corresponding to modelling the wing 
by an elastic axis with variable stiffness properties along the span 
for both bending and torsional moments . An influence coefficient con- 
cept for plate-type lifting surface might be a more general approach, 
but this was not felt warranted without a specific application in mind . 

The purpose of this study was to develop and document a simple 
approach for static aeroelastic analysis and airframe stability and 
control considerations. This documentation is a necessary prelude to 
more serious considerations on actual configurations such as new 
current Navy aircraft like the F-14 and S-3. 
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II. METHODS USED 



A. GENERAL 

In a static aeroelastic analysis there are two distinct technical 
areas requiring attention. Both the aerodynamics of the problem and 
the structural aspect must be considered. 

With the free-stream subsonic Mach number and dynamic pressure 
given, the aerodynamic loadings of the rigid wing are given as: 






(. 01 ) 



where [ A ] 



Aerodynamic influence coefficient matrix at the given 
Mach number. 




incremental loadings over the wing. 

rigid wing input angle of attack, specified at 
the control stations. 



Structural deformations on the wing due to the aerodynamic load- 
ings result in a change in the angle of attack. 

g = {cXsj (.02) 

where q = the free stream dynamic pressure, scalar constant. 

[ S 1 = the structural influence coefficient matrix. 

[a 5 ] = the angular deflections at the specified control 
points on the wing due to aerodynamic loadings . 
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Now the complete solution to the aeroelastic problem is given by: 



[lM-qLsll{j]={a r ] (.o 3a ) 

or 

ff] = N- < 3i 5 iTw — ( -° 3b) 

B. AERODYNAMICS OF THE PROBLEM 

The basic approach used in this development is similar to the 
lifting surface analysis as described by Belotserkovskii Ref. 4 . 

With the concept of a vortex lattice array of horseshoe vortices the 
wing is covered with N boxes where: 

N = NC x NS (.04) 

and NC = Number of chordwise station 

NS = Number of spanwise station 
In a manner similar to lifting line theory, as developed by Prandtl, 
Glauert, Multhopp, and Weissinger (Refs. 2 and 3), the swept bound 
horseshoe vortices were placed at the local quarter chord of each box 
and the flow tangency requirement imposed at the local three-quarter 
chord control point of each box. Figure 1 shows the geometry of the 
problem, with A and B denoting the left- and right-hand corners of 
the horsehoe vortex element and C the control point. 

As a consequence of the Helmholtz vorticity theorems (Ref. 2), 
and because inviscid flow was assumed, each horseshoe vortex ele- 
ment has constant strength, once established. The Biot-Savart Law 
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Figure 1. Horseshoe vortex lattice geometry. 




Figure 2. Finite vortex element, between A and B, 
of strength C . 
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(Ref. 2) then defines the velocity induced at a control station C (see 
fig. 2 ), for a vortex filament of finite length as: 

p 

Auj = — r (cos cc + cos^ (,05a) 

and the down wash angle as: 

= — ! — ( — ) ( COS OC 4 (05b) 

V 4rrr V ; 

If the vortex element is expanded to plus and minus infinity, the 
well known expression for downwash angle in two-dimensional flow 

results: ^ (-06) 

since cos cc = cos ^ =1.0 

As previously stated, the flow tangency requirement was applied 
at the three-quarter chord control point, and the flow was required to 
be parallel to the zero lift line at that point. Hence, 




Justification for the use of a local three-quarter chord control 
point in each vortex lattice box has not been rigorously established. 
The concept of the "rearward" aero-center at the wing three-quarter 
chord control point, from steady and unsteady two-dimensional wing 
aerodynamics, appears to be the genesis of this usage. At present 
the modelling is considered satisfactory because the theoretical load- 
ings are in accord with experimental results. It is not clear as of 
this writing why the preservation of flow tangency at local vortex- 
box three-quarter -chord control points preserves the Kutta condition 
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at the wing trailing edge, but the calculated values of local vortex 
box loadings do show the preservation of this trait. (Ref. Appendix D) 

The value of the downwash angle (-~) 3c varies at each control 

"a 

station on the wing, and there are contributions by the bound and 
trailing vortices from all of the wing stations. Hence, the flow 
tangency requirement may be expressed as a column vector as follows: 

— - <•“> 

where the matrix [a] is the aerodynamic influence-coefficient matrix 
as developed in Appendix B. 

For the vortex lattice representation with N horseshoe vortices, 
the flow tangency requirement was given by: 

[R= = N 

A 

° r [&\ = im'k] — ( - 09) 

There was a temptation to view the vortex strength at each 
control box as the pressure differential of that box. This, however, 
is not the case. Again as a result of the Helmholtz vorticity theorems 
(Ref. 2), the vortex strengths of the individual boxes may be summed 
in the chordwise direction to obtain the total vorticity at each span- 
wise station. Hence the vortex strength of a section may be "lumped" 
into one vortex equal to the sum of the individual vortices. (See Fig. 3.) 

*The development of the [ a ] matrix for spanwise stations 
only may be seen in Ref. 1. 
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Figure 3 . Contour encircling all / ' s located at the 
individual boxes . 




X 




X 



a) Field X, Y, Z b) Field X 0 ,Y, Z 

Figure 4. The effects of planform distortion. 
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From the total vortex strength at a chordal station, the section lift 

NC 

loadings can be obtained, i.e., n 2 &r- — ( - io) 

U \ 

and by application of the Kutta-Joukowski law, 

l) = 2 I* (.ii) 

a * v 

The section moment may be obtained by taking the moments about a 
reference station of the individual bound vortex elements, inasmuch 
as the bound vortex elements represent the local air loading of each 
box by the Kutta-Joukowski law. 

A two-dimensional study was made to compare the "pressure 
distribution" imposed by the vortices and the cotangent type pressure 
distribution as prescribed in theory (Ref. 2). 

Assuming that the chordwise vortex elements described a chord- 
wise pressure distribution, it appeared feasible to use a chordwise 
collocation technique for obtaining section lift coefficients and moment 
coefficients , C ? £ , given by 

c t = 2 tt /V 0 4-rtl\ v , $ 0^= 2 ( 

S lnc e = j *=!((- cose} 

9- V » c. 2 

and = -4 A_ cot § 4 A, S 

V 0 2 

it appeared that a solution for the Fourier coefficients could be found 
in the following relationship: - 



cot ^ siN2e t ••• S\m(v>-Oo i 




1 1 , • • • I 

COT SIN •• • SInGvO©^ 
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However, the existence of the A n 's assumes knowledge of the pressure 
distribution over the wing, and correct results were not obtained, 
since the vortex strengths did not represent pressure loadings. 

A check was made using two-dimensional theory to bear out 
these conclusions, and it was found that even for ten chordwise sta- 
tions, the pressure distribution was incorrect, and the aerodynamic 
center was not located at the one-quarter chord. The leading edge 
singularity was not properly represented and the first collocation 
point differed from theory significantly. This check can be seen in 
Appendix D. 

In light of these results, a straight-forward summation was used 
to obtain the section lift coefficients , i.e., 



C. COMPRESSIBILITY CORRECTIONS FOR SUBSONIC FLOW 

The compressibility corrections were obtained by the Prandtl- 
Glauert planform distortion method (Ref. 2). The applicable flow field 
equation for the velocity potential was given by 



where the subscript represents second partial derivative. By use of 




(. 12 ) 




(.13) 



the Prandtl-Glauert transformation. 




(.14) 



this reduced equation (.13) to: 



C p + 4> 4 = O 

*0*0 9VJ ^ 



(.15) 
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The effect of the Prandtl-Glauert transformation was to stretch all of 

I 

the x-coordinates by the factor 5 — . The effects can be 

M 

seen in Fig. 4. 

The net effects of the planform distortion were as follows: 

1 ) the aspect ratio was decreased; 2 ) the leading edge sweep was 
increased; 3) the taper ratio (C^/C ) was unchanged; and 4) most 
importantly, the downwash distribution { was unchanged, i.e. , 

_ d> _ 



ur 






(.16) 



Since the downwash angle distribution, , is unchanged in 

going from the distorted wing to the original wing, the aerodynamic 



loading distribution 






will be the same on both wings. This is 
an important concept in that the aerodynamic loading obtained from 
the distorted planform under incompressible flow will be the same as 
the aerodynamic loading on the original planform under compressible 
flow. Hence the aerodynamic loading is an invariant when going from 
the compressible undistorted planform to the incompressible distorted 
planform. 



D. STRUCTURAL ASPECTS OF THE PROBLEM 

Because wings are flexible or elastic, there are always changes 
in angle of attack due to the deflections of the structure. There are 
two deflection modes present, torsion and bending, both affecting the 
changes in angle of attack distribution on the wing. The effects of 
these two modes were called static aero-elastic effects and have been 
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incorporated into the program by the simple elastic-axis beam 



representation . 

The final angle of attack distribution was made up of 

two components ' anc ^ | / i« e * 

{<*<}= [<*s] + f«,} . (-17) 

Now the {o^ s *| as shown in Ref. 1, has been determined to be: 

Kl-dtsitf} 

See Appendix C for the development of the structural influence 
coefficient matrix [^s] . 



E. TOTAL SOLUTION TO THE PROBLEM 

By combining equations 0 12, , 17 , and .02, the following rela- 
tionship was formed: 

IM{^1 = 9l5l{f] 4(« t j 

“ = {«,1 ----- (•!») 

The closed form solution for aerodynamic loading distribution 
may be given as: 

= K} ( ' 03b) 

This was the final equation used in the program to obtain a chordwise 
and a spanwise loading distribution for a given surface. 
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F. CALCULATION OF STABILITY DERIVATIVES 

Five stability derivatives were investigated to check the validity 
of the program. These were: 



. , c, , k c, 



OC oC C u x p dcu 



The total wing C^ was obtained by a spanwise collocation 



technique, using the section lift coefficients. 

= Q- c ave = 2 

where =. (^) — cos © j d ^ ~ - StN G d© 
oo 

“ l=?T=EB„.».v.e 

9 V & ” 

Then \ F 1 

C, C flWe = p ) 4> B^SVViYt© S\Ni© d© 

L ave c. n _, 

oo ^ 

" = i 2 J StN StN 0 d© = jr B, 



(. 19 ) 



Since by orthogonality 
SIN Y>© 



I 



© d© 



Therefore 



=f? 

l o 



FOR Y>=t 



Tf B> 

A c. 



(. 20 ) 



AWE 

Since the assumption of linearized wing theory was made. 



C L = AC t 

oC 

whole wing , i . e . 



due to 



(-4 



C. = 



input of 1.0 rad . uniformly over the 



L «t 



The C w s were obtained by a straight-forward integration 

Oi 

technique. In effect, the section pitching moment at the quarter chord 
/v mc- 

was found to be: (^r) — ( ^1) (X- 
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Hence, total wing pitching moment is: 



C m _ = ^ i{(x-x )) 

m c/ 4 5 c <4 1 c/ A ) 

Again C.„ was evaluated by: C m =■ AC 



(. 21 ) 



v 

where AC*, is pitching moment change about the wing c /^ point 

due to R1 = 1.0 rad. input. 

With an< ^ Cy^ calculated, the was easi ly found 

by the following relationship: 



^ ^ 

CL dC L C. 

'~dL 



(. 22 ) 



From the anti- symmetric solution of the Aero Matrix (see 
Appendix B), certain rolling moment derivatives were extracted. The 
rolling moment coefficients were obtained by a collocation technique , 

similar to that used to obtain the total wing C L . 

4 I 

RM 






? = C e Cave = " i I (|) \ dl Z 



00 



*• ■= - Jr ^ 2 r\Q cos © s\^0 d 0 



O v\- \ 



= - A 2 \ S ' N ^ s ' w ^ © - \e> ^2 

® y\- 1 



where 



T\ 

I 



Therefore 



S\M Y\G SlU 



C = - 



2 eAe 



■f 



FOC W-2 

Q it Ti ^ 2. 



B-i 



7^ 



Ifo C, 



(. 23 ) 



"AME 
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The roll damping derivative, , was then obtained by placing 



the anti-symmetric angle of attack input on the wing, as shown in 

P b 

Figure 5. At the tip, the angle of attack = ^Ty = 1 rad., and 
Q P was given by: 

Kjt AC* 



C„ = 






(.24) 



The aileron control effectiveness, Cj , was found by defining 
the location of two of the boxes to coincide with the location of a 
typical set of ailerons, see Fig. 6. By inputing zero angle of attack 
everywhere on the wing except at the aileron locations and giving 
them 1.0 rad. angle of attack each, positive on one side of the wing 

and negative on the other, a rolling moment was produced. Cp was 

So. 

given by: 



C = 



^ c* 






l&r 



006 TO 






(.25) 
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Figure 5 . Angle of attack distribution for a constant 
rolling velocity P. 
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III. RESULTS AND DISCUSSION 



A. LONGITUDINAL 

The results from the symmetric portion of the program were 
checked by using longitudinal stability derivatives (i.e. lift-curve 

C /*"* d C 

L , moment-curve slope ^yy\ , and ~T ~ 7 r ~ )• To insure 
* <* 6 

that changes in leading edge sweep and Mach number are handled 
correctly in the program, checks were made on C L ^ an< ^ * 

These checks were made for a rigid wing with an aspect ratio of 6.0 
and a taper ratio of 0.3. 

Figure 7 shows the effects of leading edge sweep on lift-curve 

slope C L , and moment-curve slope C . These curves were 

* a 

obtained by inputing a "do loop" in the program and varying the lead- 
ing edge sweep -A- L£ from 0 to 40 degrees. These results agreed 



quite well with theory, i 



.i.e. (C L ^ 



(Cj. 






*■*- V i-m^c OS "a c/4 



... (.26) 



where the Mach number is held constant and the sweep varied. 

Figure 8 shows the effects of the compressibility corrections in 
the program on and . Again a "do loop" was used, 

stepping through the program for Mach numbers of 0.0 to 0.9. These 
results also agreed with theory, i.e. equation (.26), where the sweep 
is held constant and the Mach number varied. 
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oP° 




Figure 7 . Lift-curve slope plotted as a function of leading edge 
sweep at a Mach number of 0.5. 
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AR=6 0 




Figure 8. Lift-curve slope and moment-curve slope plotted as 

a function of Mach number with a leading edge sweep 
of 20.0 deg. 
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It should be noted that the variation, with increasing 



leading edge sweep, is in a direction to shift the aerodynamic center 

dc^ 



aft. Figure 9 shows this more clearly since the value of 



dC u 



about the wing corresponds to the distance, in wing m.a.c. 

lengths, that the aero center is from the reference axis. A positive 
value of 



dC 

^ indicates that the aero center is forward of the 



dC L 

reference point. 

B. LATERAL 

In order to document the anti- symmetric portion of the program, 

two lateral stability derivatives were investigated, roll damping Cj^ 

and roll due to aileron deflection Cg . The aeroelastic effects were 

So, 

checked at the same time by choosing an example from Bisplinghoff , 
Ashley, and Halfman, Ref. 3. A typical jet transport with a straight 
taper and straight elastic axis defined at 3 5%c is used. The F& for the 
wing was 6. 15 and the taper ratio was 0.444. The bending and torsional 
stiffness properties of the wing were taken from Ref. 3, Fig. 2-29; then 
they were scaled to our model by (E.l^= and reproduced 

here in Fig . 12 . 

The checks were made at a constant Mach number while varying 
the dynamic pressure (q in p.s.i.) from 0.0 (rigid wing) to 18.0, corres- 
ponding to varying the altitude. The effect of aeroelasticity upon roll 
damping for this wing, with an unswept quarter chord, is to increase the 
damping with an increase -in q. Note that the sign on roll damping is 
negative since the moment that is developed is opposite to the direction 
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AR=6.0 
TR=0.3 
M = 0.5 




Figure 9. Static margin indicator plotted as a function of leading 
edge sweep at a Mach number of 0.5. 
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of roll. Figure 11 shows the roll damping derivative plotted as a function 
of q. The results of aeroelasticity upon the rolling moment due to aileron 
deflection can be seen in Fig. 10. Note that the change in sign corres- 
ponds to the control reversal dynamic pressure, q r . The q f was found 
to be 11 . 8 p. s .1. , while it was shown to be 11.39 p . s .i. in Ref. 3 , by 
another method; this is a difference of approximately 4% . 



31 



AR=6.154 

TR=0.444 




Figure 10. Rolling moment due to aileron deflection plotted as 
a function of free stream dynamic pressure at a 
Mach number of 0.5. 
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'I 1 I I L_ 

0 4 8 12 16 



q (ps.i.) 



Figure 11. Damping in roll derivative plotted as a function of 
free stream dynamic pressure at a Mach number 
of 0 . 5 . 
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X 



Figure 12. Bending and torsion stiffness curves. 
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IV. CONCLUSIONS AND RECOMMENDATIONS 



This computer program can provide a 'limited' capability to make 
stability and control analyses of modern Navy aircraft. The program, as 
it now stands , is limited both to wings structurally compatible to the 
elastic-axis concept and to simple planforms . However, slight modifi- 
cations could easily handle arbitrarily shaped wings and even a variable 
sweep configuration, such as the F-14. 

The study is not complete at this point; a logical addition to the 
program would be to add the influence of the body radius to the Aero sub- 
routine. This would improve the value of the program. If felt warranted, 
plate theory could be used for part or all of the structures matrix. 

A follow-on study now seems appropriate to make the program more 
attractive and useful to design groups. 
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APPENDIX A. ASSUMPTIONS 



In writing this program, many theories with basic underlying 

assumptions were used. Most of them are associated with thin airfoil 

theory, or the structure of the wing. These assumptions, most of them 

listed by Gray and Schenk in Ref. 1, are as follows: 

1. Potential flow about the wing with no boundary layer effects , 
i.e. inviscid, and without compression shocks, i.e. below 
transonic mach numbers . 

2. The thickness ratio of the wing is small. 

3. The lift-curve-slope, CL = 2TC ) i.e., two-dimensional 

4. Changes in camber due to bending and torsion of the wing 
are neglected. 

5. The changes in angle of attack due to structural deformations 
are uniform along each spanwise station. 

6. The elastic deformation of the control surfaces is the same 
as that of the adjoining wing. 

7. The angles of structural deformation are small such that 

SIN © = TA.U © = © j COS © = \.0 

8. Drag-load effects to deformation of the wing were neglected. 
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APPENDIX B. DERIVATION OF THE AERO SUBROUTINE 



The basic development of the downwash matrix or aerodynamic 
influence matrix was done by Gray and Schenk in Ref. I. An expansion 
was made here to include swept-bound vortices and a chordwise, as well 
as spanwise, distribution of horseshoe vortices. 

The lift or aero loading distribution over an airfoil can be charac- 
terized by a system of horseshoe vortices, each of these vortices con- 
tributing to the upwash or downwash at points all over the wing. Figure 
B1 shows two horseshoe vortices combined to make a wing. Conventional 
right hand rule indicates the sense of circulation associated with each 
horseshoe vortex. Quite obviously, the more horseshoe vortices used to 
approximate a wing, the more accurate the lift distribution obtained. 

With the wing divided up into NC X NS stations, the horseshoe 
vortices were placed at each quarter chord of the box, and the horseshoe 
vortex lattice, as shown in Fig. 1, was arrived at. Now the boundary 
conditions or flow tangency requirements , for simplicity of calculations , 
were imposed at the three-quarter chord * for each station or box. The 
boundary condition implies that the downwash angle ( ^) at each con- 
trol station is equal to the geometric angle of attack. The downwash vel- 
ocity U7 at each station is due to the sum of the downwash velocities in- 
duced by each horseshoe vortex over the entire wing. 

*The chord for each box is equal to the local wing chord divided 
by NC. 
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1 




Figure Bl. Simplified horseshoe vortices system for 
wing lift distribution. 





X 

Figure B2 . Geometry for downwash velocity calculation. 
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Given one horseshoe vortex and one control point on the same 



wing, as shown in Fig. B2 , the downwash velocity Aur was calculated 
as follows: 

R* = V(y c -y A f + (x c -y a ) z •> R 6 =V ( y- + ( v c -x B y 

SmQ = (x c -x a Vr a s\nc^ = (x c -x b Vr b 

cos © = Cy c - y*) /Ra cos $ = ( y c - y B V 

Ca B = V (y a -y A f + (x B -y A ^ 

5\N (Xg-X/jJ)/ c OS / J^AB 

OC= G - X' 4 ~ ^ 

S\N C< =r SIN Q COS X — COS© SIN X 

COSOC = COS© CO S ^ + S\N G SIN ^ 

cos ^ = coscj) COSjf 4- S\N § S\N X 



The left hand vortex contribution to the downwash velocity was given by: 

(B.01) 



A , , _ r o + sin q) 

AW,— ~A — 7 , — 

An (y c -yJ 



The right hand vortex contribution was given by: 

A , -r ( i-v sin^ 

Au) = ^ — 7 r~ 

R 4 n (y c -ij B ) 

The bound vortex contribution was given by: 

„ . , r (cos#- COS0) 

AOJ r =l 



4r< 



£ A sikioC 



(B. 02) 



(B. 03) 
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If the control point happened to be on a straight line extension of 
the bound vortex, both numerator and denominator would go to zero 
simultaneously . The contribution from the bound vortex would be zero 
and a check was placed into the program and if this situation were to 
develop, AuJ^ would be set to zero. 

The total contribution to the downwash velocity due to the one 
horseshoe vortex on the same wing was: 

Au> = ku>44 Auj u + 

or &uJ= ^ ^ (costt- cosfe’)! (b. 04 ) 

4rtUy c -%) 1 

The downwash velocity contribution at control points all over the 
wing, due to a horseshoe vortex system, was calculated in a similar 
manner. These can be thought of in matrix format as follows: 




A,U> u + 

4 AW w 



4 

4 

i 

l 

i 

i 

4 



4 AU3 IN 



. 4 &00 






• 4 Aw 



NN 



where the subscripts indicate the downwash at station i due to a horse- 
shoe vortex at station j, etc. Now by factoring a ]~* from each term 
in the proper manner: 



W=zUlW 



then the following relationship follows: 
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(B. 05) 



tsl- 



(B.06) 



Now 



{vi=K] > = — <“-» 7 > 



The A matrix is the aero-influence matrix as used in the program. 
The computations are carried out for contributions from each wing, i.e. 



u«i i ua 

For symmetric loads: l \ 1 

For anti- symmetric loads: [&M=U r H a /I 

ah 



(B.03a,b) 



Then I Of /l = t AS or* A&"] { ^r-\ 

^ 1 } > K.v*. 



S\oe 



Also contained in the Aero subroutine is the Prandtl-Glauert plan- 
form distortion correction for subsonic compressibility. In essence, all 

Z\' /z 

of the x distances are stretched by a factor of * , defined as 

. Therefore throughout the subroutine replaces . 

(i.e. instead of C X c - X ^ ). 

The 'lifting-surface' representation of a wing has been written up 
in FORTRAN IV for use with the IBM 360-67 computer. The program was 
designed to handle planforms of low to high aspect ratio, swept leading 
edges, and a linear taper ratio. The choice of the number of stations on 
the wing, chordwise and spanwise, was left arbitrary for easy compati- 
bility with most wings and for almost any flap and aileron arrangement. 
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So as to make the program more flexible, subroutines were used 
for each phase of the calculations. Subroutine Geom is used for locating 
the horseshoe vortices and their control points; subroutine Aero is used 
for obtaining the aerodynamic influence (or downwash) matrix; and sub- 
routine Struct is used for obtaining the structural influence coefficient 
matrix. With the subroutine packages set up as they were, changes to 
the program can be easily handled simply by making small changes in the 
input or the main program. 

Sample calculations were made for arbitrary wing geometries . The 
results of these calculations can be seen in the main text. 

The desired results were obtained in these calculations. The 
corrections for subsonic compressibility, sweepback, and static aero- 
elastic effects have been checked and do give the desired accuracy. 



42 



APPENDIX C. DERIVATION OF STRUCTURES SUBROUTINE 



As in Appendix B, the [.S] , structures matrix or elasticity matrix, 
was developed by Gray and Schenk in Ref. 1. An expansion to fit 
'lifting-surface' theory vice 'lifting-line' theory was made for this 
program . 

The problem was basically stated as: Given a loading over a wing, 
find the deformation imparted to the wing and apparent change in angle 
of attack. From the previous subroutine, 



The following relationship was obtained from this subroutine: 

fVsi = <3- (-02) 

The assumption has been made that the wing acts like a beam with 
the elastic axis representation. (See Refs. 1 and 4.) The beam bending 
and torsional moments have been resolved about the elastic axis at each 
spanwise station of the wing. Although the elastic axis may be arbitrary, 
it was thought of as a straight line for simplification here (see Fig. Cl). 

The wing is divided evenly in the spanwise direction with each of 
the chordal strips having width 2h. Therefore for the wing with an 

assumed unit value of semispan, h = 1/2 (NS) (C.01) 

Positive direction for torsion and bending moments are shown in Fig. Cl 
relative to the elastic axis orientation. 
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Figure Cl. Elastic axis representation of the wing with 

the torques and moments at the station elastic 
axis points . 
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for the 



The span loads per unit span length are defined as • 

-r 

N vortex lattice boxes in a numbering sequence starting from N = 1 at 
the most forward, inboard box proceeding aft to the trailing edge at 
N = NC, then stepping to the next outboard set of boxes at the leading 
edge box for N = NC+1. Continue in this manner until the Nth box at 
the trailing edge of the tip is reached. 

The individual lift force contribution towards the bending of the kth 
wing station may be expressed as: 

^U2-S.*) qln (C.02) 

The presence of the Kroecker delta takes into account that only one-half 
of the air loading at the kth station contributes to the bending (and torsion) 
moment at the kth station. 

Therefore the bending moment at the kth station may be expressed 
as: 

NS 

M. = (.%) + 2 ] ALj (c . 03) 

or in matrix notation as: 

= % \.MA\ (C.04) 



where 


= tvs* \ 


column vector 




{$1 - X * ' 


column vector 




= MS * N 


rectangular matrix, and 
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[ami = e? 



'X 

o 

o 

o 



\a 


8 1 12 I l £ 1 - 


h/~ 


|. 1 _L _i 

, 4. 1 A 1 t o \ 


h-N 

° j 


1— - 2 I 

l /n \ A a 


i 'Z J 2 ' L 

fo_l Li 4- — — • 

I 1 l | \L | 


i < 1 


i i 1 
1 1 1 


• 1 ' -i— U- 

1 : 1 : 1 ! i J 



The AM matrix was partitioned into NS x NS elements where each ele- 
ment is a (1 x NC) row matrix made up of constants as indicated, i.e. 

J l /z 1 = 1_ '4 '/z • • • ■ 

~l 7 



The torsion moment at the kth station may be expressed as: 



T*= SIAT]^} 



(.05) 



where 



[*T ] = * 



AX, 



AX V 



_o_ 

O 



WC |2 AX^ KC4% * * * 2 Ax, ?Kc j 
o I 



AX 



“T 



I 



T 



and 



where 



^ X c/ A )j 



^ I ^ J ^ Ni ~ NIC * NS 
The zero partition elements to the left of the lead diagonal correspond 

to inboard loadings that do not contribute to the moments at stations 
outboard . 
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The torsional and bending moments as developed here were perpen- 
dicular and parallel respectively to the streamwise flow at the elastic 
axis reference points. In order to reorient them relative to the elastic 
axis , they were put through the following rotation matrix: 



( 




cosAe* 




M 


S\N A 6a . 


cosA 


M 



(C . 06) 



Now the moment and torque distributions have been resolved to be the 



following relationships: - 



The change in angle of attack in the streamwise direction was ob- 



tained by the relationship (see Ref. 1): 






.= 5 



vn 



Mds 



El 



+ 



t Ids 

GJ 



(C . 08) 



o w o 

where m equals the beam bending moment per unit pitching moment at 
the station , and t equals the torsion about the elastic axis per 

unit pitching moment. ds => r __ (C.09) 

Co s ^ 

The above integral equation was changed to an algebraic equa- 
tion by summing the contributions at each station. 



(C . 10) 





^ 2 Yr»2, ^ 


(E3\cosJL. 


(n\cosA t 


"tu Tt ?>i 




CGJ),cosJ\., 


(Gj^cosA^ 
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The influence coefficients m^ and tjj were easily defined by the geom- 
etry of the problem and found to be as follows: 



and 



t ij = 



^-sin A • where the point i is outboard of the point j 
0 where the point i is at or inboard of the point j 

( cos for station i outboard of station j 

0 for station i at or inboard of station j 



In matrix form the m^'s and t^'s were written as follows: 



H= 



M- 



2 


- simA 1 


"SIN 


0 


- 

2 


- SIN A-*, 


0 


0 


_ s\u As 


1 


1 

1 


Z 

1 


cosAi 

2 


cosA-z. 


i 

cos Jl z 


O 


cosAi 

2 


COsA-j 

cosAs 


O 


O 


2 



Now 



{“4 = + [■ 4[o&a]{t '} 



2 * 
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Thus the desired form was achieved and the S matrix was: 



[s]- jjco s A][AM] - [s\n a] [at]] 

+ Ulf^j][[s^A][AV\] + [co5A][AT]]] (C . 12) 



The final resultant matrix equation was: 



The computer program is normalized such that the wing semispan 

has unit length. However / an actual wing will have a finite wing span 

of say b/2 inches, will operate at a dynamic pressure of q p.s.i., and 

2 

will have stiffness properties of El and GJ in units of lb-in . 

Recognize that in a scaling analysis the structural twist of the 
model will equal the structural twist of the prototype and that they will 
operate at the same dynamic pressure. However, the aerodynamic load- 
ings will differ by a factor of \ 



where 







i.e. 



hl= i {% 



Therefore by equation (.02) we can scale 



(C . 13) 
(C . 14) 



1^]^= "Alslp ^ 

but analysis of equation (C.12) reveals that Is] OC WEI , V6J 

~ -a. ^ ) 



[s]„ . 



U1 P a 



but 









(£V)p= 7t A (Gl] m 
(G = A a (GJU 



- (C.15) 
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APPENDIX D. TWO-DIMENSIONAL WING PRESSURE 



LOADING CONVERGENCE CHECK 



In this appendix, a check was made to ascertain the validity of 
viewing the vortex strengths of the control boxes as pressure differentials. 
A comparison was made for the C^'s obtained by this assumption, using 
a collocation technique for obtaining a chord wise pressure distribution, 
and the theoretical Op’s with a cotangent type pressure distribution. 

A two-dimensional section with unit chord was chosen for this 
check. The section was divided into N equal chordwise stations with 
the vortices placed at the quarter-chord of each station and the control 
points at the three-quarter chord of each station (see Fig. Dl). 

Since C > Q f Q represent the small vortex strengths concentrated 
at the quarter chord for each station, the total vortex strength was given 
by: 



Now 






i a Pl= M = £l£ 

Y 1 AC C/N 



and 



* c r,= -£^7i 



2 m 



(D.01) 



'/ 2 pv z 0/V> v V 

The principle of an infinite vortex (Ref. 2) was used for the down- 

w _ _j t r\ 

wash angles, i.e. \j ~ Z n R '•V ' or 




(D. 02) 



where — 

V 



OC at the three-quarter chord control point, and 







(D . 03) 
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V 




Figure Dl. Geometry for two-dimensional study of 
chordwise pressure distribution. 
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, the local quarter chord point, 



with 




A-’ - 


a 75 


and 


*c- = 


A^- 


0.25 

M 



, the local three-quarter chord point. 



Hence {&Cp| = 4 7\ N t Bl {<*} (D.04) 

with OC =1.0 radian for each station. 

Now the theoretical & for two dimensions was given by 

(Ref. 2): AC- = 4 

? M0 

using the coordinate transformation: 

f = -i ( l - Cos <$>) 

Q = cos' [ l- 2(|M 

Since the check was made at OC = 1.0 radian and a unit chord, the 

obtained will be equivalent to and should be approximately 

C m 

27\ . The aerodynamic center given by: AC = — (D.07) 

should be at 0.25c or just 0.25 in this case. 

The tabulation in Fig. D2 clearly shows that the correct results 
were not obtained by the method as described herein. 



(D.05) 



(D . 06) 
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Figure D2. A tabulation of the two-dimensional wing convergence check. 
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SUBROUTINE ELU TR I-D I AGONAL I 2ES THE INPUT SQUARE MATRIX FOR SOLUTION 
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420 AM { 1 1 J ) =0 • 0 

CG 4 JO 1=1 ,I\!Z 
DO 430 J=1,NZ 
DO 4 30 K=1 f NS 

B ( I , J) =MI ( i ,K) *V (K, J >+3(I, J > 

4 30 AM( I , J) = TI C I tK)*U(K v J) + AMU » J) 
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